Cardiac mechanics Benchmark - Problem 3

This code implements problem 3 in the Cardiac Mechanic Benchmark paper

Land S, Gurev V, Arens S, Augustin CM, Baron L, Blake R, Bradley C, Castro S, Crozier A, Favino M, Fastl TE. Verification of cardiac mechanics software: benchmark problems and solutions for testing active and passive material behaviour. Proc. R. Soc. A. 2015 Dec 8;471(2184):20150641.

import dolfin

try:
    from dolfin_adjoint import Constant, DirichletBC, Mesh, interpolate
except ImportError:
    from dolfin import DirichletBC, Constant, interpolate, Mesh

import pulse

geometry = pulse.HeartGeometry.from_file(pulse.mesh_paths["benchmark"])

2021-05-13 08:28:25,025 [72707] INFO     pulse.utils:
Load mesh from h5

# Create the material
material_parameters = pulse.Guccione.default_parameters()
material_parameters["CC"] = 2.0
material_parameters["bf"] = 8.0
material_parameters["bfs"] = 4.0
material_parameters["bt"] = 2.0

activation = Constant(0.0)
material = pulse.Guccione(
    parameters=material_parameters, active_model="active_stress", activation=activation
)

# Define Dirichlet boundary. Fix the base_spring
def dirichlet_bc(W):
    V = W if W.sub(0).num_sub_spaces() == 0 else W.sub(0)
    return DirichletBC(
        V, Constant((0.0, 0.0, 0.0)), geometry.ffun, geometry.markers["BASE"][0]
    )

# Traction at the bottom of the beam
lvp = Constant(0.0)
neumann_bc = pulse.NeumannBC(traction=lvp, marker=geometry.markers["ENDO"][0])

# Collect Boundary Conditions
bcs = pulse.BoundaryConditions(dirichlet=(dirichlet_bc,), neumann=(neumann_bc,))

# Create problem
problem = pulse.MechanicsProblem(geometry, material, bcs)

# Solve problem
pulse.iterate.iterate(problem, lvp, 15.0, initial_number_of_steps=50)
pulse.iterate.iterate(problem, activation, 60.0, initial_number_of_steps=50)

2021-05-13 08:28:25,420 [72707] INFO     pulse.utils: Iterating....
2021-05-13 08:28:25,421 [72707] INFO     pulse.utils: Current control: 0.000
2021-05-13 08:28:25,421 [72707] INFO     pulse.utils: Target: 15.000
2021-05-13 08:28:47,766 [72707] INFO     pulse.utils: Iterating....
2021-05-13 08:28:47,767 [72707] INFO     pulse.utils: Current control: 0.000
2021-05-13 08:28:47,767 [72707] INFO     pulse.utils: Target: 60.000

([Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 294),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 320),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 344),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 370),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 398),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 426),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 454),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 615),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 645),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 804),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 836)],
 [Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 292),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 318),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 342),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 368),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 396),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 424),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 452),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 613),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 643),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 802),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 834)])

# Get displacement and hydrostatic pressure
u, p = problem.state.split(deepcopy=True)

endo_apex_marker = geometry.markers["ENDOPT"][0]
endo_apex_idx = geometry.vfun.array().tolist().index(endo_apex_marker)
endo_apex = geometry.mesh.coordinates()[endo_apex_idx, :]
endo_apex_pos = endo_apex + u(endo_apex)

print(
    ("\n\nGet longitudinal position of endocardial apex: {:.4f} mm" "").format(
        endo_apex_pos[0]
    )
)

Get longitudinal position of endocardial apex: 11.8550 mm

epi_apex_marker = geometry.markers["EPIPT"][0]
epi_apex_idx = geometry.vfun.array().tolist().index(epi_apex_marker)
epi_apex = geometry.mesh.coordinates()[epi_apex_idx, :]
epi_apex_pos = epi_apex + u(epi_apex)

print(
    ("\n\nGet longitudinal position of epicardial apex: {:.4f} mm" "").format(
        epi_apex_pos[0]
    )
)

Get longitudinal position of epicardial apex: 15.4904 mm

V = dolfin.VectorFunctionSpace(geometry.mesh, "CG", 1)
u_int = interpolate(u, V)
mesh = Mesh(geometry.mesh)
dolfin.ALE.move(mesh, u_int)

from fenics_plotly import plot, set_renderer
set_renderer('notebook')
fig = plot(geometry.mesh, color="red", show=False)
fig.add_plot(plot(mesh, opacity=0.3, show=False))
fig.show()